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Abstract. We present the resuhs of Monte Carlo simulation for a Kondo lattice model in 
^ , which itinerant electrons interact with Ising spins with spin-ice type easy-axis anisotropy on 

$_^ ' a pyrochlore lattice. We demonstrate the efficiency of the truncated polynomial expansion 

jy^ ' algorithm, which enables a large scale simulation, in comparison with a conventional algorithm 

using the exact diagonalization. Computing the sublattice magnetization, we show the 
C^ ' convergence of the data with increasing the number of polynomials and truncation distance. 

a 

jj ■ 1. Introduction 

Interplay between localized spins and itinerant electrons has been one of the major topics in the 
field of strongly correlated electrons. It triggers various interesting phenomena, for instance, a 
^ ■ variety of magnetic orderings induced by electron-mediated spin interactions, such as Ruderman- 

'nJ" ', Kittel-Kasuya-Yosida (RKKY) interaction [1] and the double-exchange interaction [2]. The 

impact of the spin-charge interplay is not limited to the magnetic properties, but also brings 
about peculiar electronic states and transport phenomena, such as non Fermi liquid behavior 
in the quantum critical region in rare-earth systems [3] and the colossal magneto-resistance in 

$^ ' perovskite manganese oxides [1]. 

Recent studies on metallic pyrochlore oxides have opened yet another aspect in the field, 
namely, the geometrical frustration. The effect of geometrical frustration on the interplay 
between itinerant electrons and localized spins has attracted much interest, and extensive 

k> I number of studies on these compounds have been reported. For example, interesting features 

^ ' were reported in Molybdenum compounds, such as an unconventional anomalous Hall effect [5] 

and emergence of a peculiar diffusive metallic phase [B]. A characteristic resistivity minimum 
was also observed in an Iridium compound [?]• For these phenomena, the importance of local 
spin correlations inherent to the strong frustration has been suggested, but comprehensive 
understanding is not reached yet. 

One of the authors and his collaborator recently reported unbiased Monte Carlo (MC) 
calculations of a Kondo lattice model on a pyrochlore lattice [H O [TO]. In their studies, 
however, the accessible system size N was limited to small sizes because the MC simulation 
was performed by a conventional algorithm using the exact diagonalization (ED) [11], in which 
the computational amount increases with 0{N'^). Larger size calculations are highly desired to 
further discuss the peculiar magnetic and transport phenomena in the frustrated spin-charge 
coupled systems. For this purpose, here we apply another faster algorithm, the polynomial 



expansion method (PEM) |121 I13j . We test the efficiency of the algorithm for a variant of 
Kondo lattice models on a pyrochlore lattice. 

2. Model and method 

We here consider a Kondo lattice model on a pyrochlore lattice, whose Hamiltonian is given by 

^ = -i E (cUci,<x + H.c.)-J^S,-^,. (1) 

Here, Cj^o- (c^^) denotes an annihilation (creation) operator of electrons at site i with spin 
o"(=t, ^|.); Sj and <tj represent the localized spin and itinerant electron spin, respectively. The 
model is defined on a pyrochlore lattice, three-dimensional frustrated lattice consisting of a 
corner-sharing network of tetrahedra. The sum {i,j) is taken over the nearest-neighbor (n.n.) 
sites on the pyrochlore lattice. We set i = 1 as the energy unit. 

We assume that the localized spins Sj are Ising spins with |Sj| = 1, whose anisotropy 
axes depend on the four-sublattice sites on each tetrahedron: the axes are set along the (111) 
directions, namely, the directions connecting the centers of neighboring tetrahedra. The situation 
is the same as in the spin ice |141ll5j. Although there is no bare interaction between the localized 
spins in the model ([1]), the spins communicate with each other through an effective interaction 
mediated by the kinetic motion of electrons (RKKY interaction). Similar to the spin ice case, 
when the n.n. interaction is dominantly ferromagnetic (FM), a local spin configuration with two 
spins pointing in and the other two pointing out (two-in two-out) is favored in each tetrahedron. 
On the other hand, when the interaction is dominantly antiferromagnetic (AFM), all-in or all-out 
configurations become energetically stable. We note that the sign of the coupling J is irrelevant 
since the Hamiltonian is unchanged for J — t- —J and S, — ?■ — Sj. 

We apply MC calculations with PEM to the model ([1]). In the PEM, the density of states 
(DOS) for electrons under a spin configuration is obtained using the Chebyshev polynomial 
expansion |12j . which then is used to calculate the action in the weight for the given spin 
configuration. This algorithm reduces the computational cost compared to a conventional MC 
method based on ED [11] . We also implement the truncation method which further reduces 
the computational cost [13]. We carry out the truncation by a real space distance, not by a 
magnitude of the matrix element in the original scheme; namely, we introduce a truncation 
distance defined by the Manhattan distance from a flipped spin in calculating the Chebyshev 
moments. The total computational amount for each MC step is largely reduced from 0{N'^) for 
the conventional algorithm to 0{N) [13]. This enables us to access to larger system sizes. 

The calculations are done with varying the number of polynomials 15 < m < 50 and 
truncation distance 3 < d < 7. In the following, we take J = 2 and restrict ourselves to 
two typical electron densities; a low electron density rig = Yli ai'^ia'^i,^) /-^ ~ ^.03 (the chemical 
potential is fixed at // = —5.9), for which the n.n. RKKY interaction is FM, and an intermediate 
electron density rif, ~ 0.35 (/i = —1.3), for which the n.n. RKKY interaction is AFM. All the 
following calculations were done for A^ = 4 x 4^ with periodic boundary conditions, typically 
with 2700 MC steps after 500 steps of thermalization. The results are divided into three bins 
of 900 steps to estimate the statistical error with 100 steps intervals in between each bin. MC 
simulation costs about 15 hours with 8 cpu parallelization of Intel Xeon E5502 processors [16]. 

3. Results 

Here, we show the MC results for the square of the sublattice magnetization, m^, with 
different m and d. We calculate ml by the diagonal component of the spin structure factor 
S"/3(q) ^ En,z(SnSf ) exp(iq • rni)/N, as ml = 4S""(q = 0)/A. Here, n and / are the indices 
of tetrahedra and a, /3 = 1, 2, 3, 4 are the indices of sites in a tetrahedron. 
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Figure 1. MC results 

by the PEM for ml at an 
intermediate density Ue ~ 
0.35. We take d = 6 in (a) and 
m = 40 in (b). The results 
and errors of the ED method 
are shown by horizontal solid 
lines and shades. 
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Figure 2. MC results by 
the PEM for m,^ at a low 
density rif. ~ 0.03. We take 
d = 6 in (a) and tti = 40 in 
(b). The results and errors 
of the ED method are shown 
by horizontal solid lines and 
shades. 



Figure [T] shows the PEM results of m}g at an intermediate density rig ~ 0.35 for three 
typical temper atures(T). Figure [T^a) shows m dependence at d = 6 and figure mb) is for d 
dependence at m = 40. At this electron density, the localized spins align in the all-in/all-out 
configuration (alternative arrangement of all-in and all-out tetrahedra) at low T since the n.n. 
RKKY interaction is dominantly AFM. The horizontal solid lines denote the results obtained by 
the ED MC method. The ED result at T = 0.07 shows rra^ > 0.8, which implies the system to be 
in the AFM ordered phase. At T = 0.09, m^ ~ 0.1, which indicates a disordered paramagnetic 
state. The intermediate T = 0.08 is presumed to be around the critical point. In figure [H^a), 
the PEM results at T = 0.07 and T = 0.09 converge to the ED results when m > 30. The data 
at r = 0.08, however, show slower convergence: It appears to require m > 35 ~ 40 for reliable 
calculation. The convergence as to d shows a similar tendency, as shown in figure [T]^b); d > 4 is 
sufficient at T = 0.07 and T = 0.09, while a larger d> 5 appears to be necessary at T = 0.08. 

In lower density regions, the convergence becomes poorer because of several reasons. One is 
simply because the relevant T range becomes lower; the RKKY energy scale gets smaller for lower 
density. Another reason is that the Fermi energy comes close to the band edge (bottom); the 
PEM is an expansion technique of DOS, and hence, larger m is necessary for good convergence 
when DOS changes rapidly as in the band edge. Furthermore, when the electron density is small, 
a fluctuation of the density in MC measurements will harm the precision; in our calculations, 
ne has typically an error of An^ ~ 0.01, which might affect the results when Ue ~ Aue- 

Such situation is illustrated in figure [2] for an extremely low density rig ~ 0.03. In this region, 
the lowest T state is characterized by a EM order of two-in two-out tetrahedra with aligning 
the net moment of each tetrahedra along a (100) direction. Note that the T range in figure [2] 
is much smaller than in figure [TJ The data show much poor convergence to the ED results: 
Nevertheless, m > 40 and d>l will give converged results even in this extreme case. 



Our result indicates that reasonable convergence to the ED result is achieved by taking 
m ~ 35 — 40 and d ~ 5 — 7 in the PEM in a wide range of the electron density and temperature. 
The convergence is considerably slower than in the previous study for a double-exchange model 
on the three-dimensional cubic lattice |17j . which showed sufficient convergence within m < 8. 
This is partly due to the difference of the magnitude of J. In the previous study, J was taken 
to be infinity, which greatly simplifies the band structure, while we take a much smaller value 
J = 2 in the present study. The resultant complicated band structure requires a larger number 
of m to reproduce its fine details. Another possibility is the effect of geometrical frustration 
in the present model. In general, the frustration leads to degeneracy among different spin 
configurations in a small energy window. For such situation, a small systematic error in the 
calculation of the MC weight could lead to considerable errors in observables. Furthermore, 
the frustration suppresses all the energy scales, which might also require much larger efforts to 
obtain converged data. 

4. Summary 

To summarize, we examined the efficiency of the Monte Carlo simulation using the polynomial 
expansion method for a frustrated Kondo lattice model with spin-ice type Ising spins. By 
comparison with the results by the exact diagonalization method, our result indicated that the 
polynomial-expansion Monte Carlo calculation with m ~ 35 — 40 and d ~ 5 — 6 gives consistent 
results in a wide range of electron density and temperature. Larger m and d are necessary 
for convergence near the critical temperature and for very low electron density. Our results 
demonstrate that the polynomial expansion method is a practical method for investigating the 
physics of spin-charge coupled systems even when J is much smaller than the bandwidth. Monte 
Carlo study of the phase diagram by systematic analysis up to larger system sizes is in progress. 
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